function example2

%  Solves the BVP:
%       y'' + p(x)y' + q(x)y= f(x)   for xL < x < xr  '
%  where
%      y(xl) = yL  and y(xR) = yR

%  p=-x^2/ep, q=-1/ep, f=0
%  xL=0, yL=1  and  xR=1, yR=1

% clear all previous variables and plots
clear *
clf

% set boundary conditions and parameters
	ep=0.01
	xL=0; yL=1;
	xR=1; yR=1;

% loop used to increase N value
for in=1:3

	% set number of points along the x-axis
	if in==1 
		n=10
	elseif in==2
		n=20
	elseif in==3
		n=120
	end;

	% generate the points along the x-axis, x(1)=xL and x(n+2)=xR
	x=linspace(xL,xR,n+2);
	h=x(2)-x(1);

	% calculate the coefficients of finite difference equation
	a=zeros(1,n); b=zeros(1,n); c=zeros(1,n); z=zeros(1,n);
	for i=1:n
		a(i)=-2+h*h*q(x(i+1),ep);
		b(i)=1-0.5*h*p(x(i+1),ep);
		c(i)=1+0.5*h*p(x(i+1),ep);
		z(i)=h*h*f(x(i+1));
	end;
	z(1)=z(1)-yL*b(1);
	z(n)=z(n)-yR*c(n);

	% solve the tri-diagonal matrix problem
	y=tri(a,b,c,z);
	y=[yL, y, yR];

	% plot the solution
	if in==1
		plot(x,y,'-.r','LineWidth',1)
		hold on
		legend(' N = 10',3);
		
		% define title and axes used in plot
		box on
		xlabel('x-axis','FontSize',14,'FontWeight','bold')
		ylabel('Solution','FontSize',14,'FontWeight','bold')
		title(['BVP: Example 2 with \epsilon = ', num2str(ep)],'Color','k','FontSize',14,'FontWeight','bold');
		% Set the fontsize to 14 for the plot
		set(gca,'FontSize',14);
		set(findobj(gcf,'tag','legend'),'FontSize',14,'FontWeight','bold');		
		
		pause
	elseif in==2
		plot(x,y,'--b','LineWidth',1)
		legend(' N = 10',' N = 20',3);
		set(findobj(gcf,'tag','legend'),'FontSize',14,'FontWeight','bold');
		pause
	elseif in==3
		plot(x,y,'--k','LineWidth',1)
		legend(' N = 10',' N = 20',' N = 120',3);
		set(findobj(gcf,'tag','legend'),'FontSize',14,'FontWeight','bold');
	end;
end
hold off


function g=q(x,ep)
g=-1/ep;

function g=p(x,ep)
g=-x^2/ep;

function g=f(x)
g=0;


% tridiagonal solver
function y = tri( a, b, c, f )
N = length(f);
v = zeros(1,N);   
y = v;
w = a(1);
y(1) = f(1)/w;
for i=2:N
    v(i-1) = c(i-1)/w;
    w = a(i) - b(i)*v(i-1);
    y(i) = ( f(i) - b(i)*y(i-1) )/w;
end
for j=N-1:-1:1
   y(j) = y(j) - v(j)*y(j+1);
end